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ABSTRACT 

We compute the molecular line emission of massive protostellar disks by solving the equation of 
radiative transfer through the cores and disks produced by the recent radiation-hydrodynamic sim- 
ulations of Krumholz, Klein, & McKee. We find that in several representative lines the disks show 
brightness temperatures of hundreds of Kelvin over velocity channels ~ 10 km s~^ wide, extending 
over regions hundreds of AU in size. We process the computed intensities to model the performance 
of next-generation radio and submillimeter telescopes. Our calculations show that observations us- 
ing facilities such as the EVLA and ALMA should be able to detect massive protostellar disks and 
measure their rotation curves, at least in the nearest massive star-forming regions. They should also 
detect significant sub-structure and non-axisymmetry in the disks, and in some cases may be able to 
detect star-disk velocity offsets of a few km s~^, both of which are the result of strong gravitational 
instability in massive disks. We use our simulations to explore the strengths and weaknesses of differ- 
ent observational techniques, and we also discuss how observations of massive protostellar disks may 
be used to distinguish between alternative models of massive star formation. 

Subject headings: accretion, accretion disks — ISM: clouds — methods: numerical — radiative transfer 
— stars: formation 



1. INTRODUCTION 

Determining the nature of the disks or other struc- 
tures that feed gas into massive protostars is one of 
the major observational challenges in radio and submil- 
limeter astronomy. Observations to date have revealed 
large-scale toroids in sub-Keplerian rotation, and in a 
few cases evidence for Kepleria,n rotation (see recen t re- 
views bvlCesaroni et al.l l2006l l2007l iBeuthej [20061 . and 
■ iBeuther et al.ll2007D . However, even the most sensitive 
' observations thus far have only been able to probe length 
scales of thousands to tens of thousands of AU from the 
embedded massive star. True accretion disks arc likely 
to begin at smaller distances from their host stars, so 
determining the nature of massive stellar disks and their 
role in the star formation process will require yet higher 
resolution. 

Observationally, this presents a significant challenge. 
Even the closest regions of massive star formation are 
~ 0.5 kpc away, so resolving structures hundreds of AU 
in size requires tenth-of-an-arcsecond angular resolution. 
A more typical distance for massive protostellar cores is 
1 — 2 kpc, requiring even higher resolution. Since the 
strong ionizing flux from a massive star rapidly alters 
or destroys its natal environment, it is only possible to 
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observe the structures from which a massive protostar 
forms while it is still embedded and accreting. Such 
embedded stars are generally obscured by hundreds of 
magnitudes of dust extinction at visual wavelengths, and 
can therefore only be observed at radio or submillime- 
ter wavelengths. Furthermore, fully understanding the 
behavior of these structures requires kinematic informa- 
tion, which can only be obtained unambiguously from 
line emission. Probing molecular lines requires signifi- 
cantly higher sensitivity than continuum observations. 

Unfortunately, no existing telescopes are capable of de- 
tecting molecular line emission at brightness tempera- 
tures hundreds of Kelvin or less with resolutions ~ 0."1. 
However, this combination of sensitivity and resolution 
is within reach of next-generation observatories such as 
the EVLA and ALMA, and studying the environments 
in which massive protostars form is one of the major sci- 
ence goals of these facilities. Fulfilling this goal will re- 
quire significant theoretical input, both to suggest what 
observing strategies are likely to be successful, and to 
make predictions which, when compared with observa- 
tions, can be used to help decide between alternative 
models. 

Recently, iKrumholz et all (|2007bl hereafter KKM07) 
simulated the early stages of collapse and fragmenta- 
tion of massive protostellar cores. These simulations 
are ideal for observational comparisons for several rea- 
sons. First, they begin from realistic, turbulent initial 
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conditions, and thus should have levels of structure and 
complexity comparable to real massive star formation 
regions. In comparison, many other simulations of mas- 
sive star formation begin fro m far simpler initial condi- 
tions with no t urbulence (e.g. lYorke fc So nnhalterl l2002t 
iBaneriee fc Pudritz 2 007? ). While these simulations are 
an excellent tool for isolating and exploring some of the 
physical problems of massive star formation, their sim- 
plicity makes them less suited for detailed observational 
comparison than more complex models. 

Second, KKM07 include radiative transfer as well 
as hydrodynamics and gravity in their simulations, in- 
cluding radiative heating by the forming protostar(s). 
As a result, the simulations make realistic predictions 
for the distribution of temperatures, as well as den- 
sities and velocities, inside a star-forming core. In 
contrast, models that neglect radiation and treat the 
gas using only a modified e quatio n of state (e.g. 
iBonnell et 311120041 : iDobbs et al.l 12005?) . or that include 
radiation only via cooling cu rves (e.g. IBaneriee et al.l 
l2006t IBaneriee fc Pudritdl2007D . do not produce reahstic 
temperature distributio ns because they n eglect heating 
by accreting protostars (|Krumholdl2006bl KKM07). 

In this paper we use the simulations of KKM07 to pre- 
dict the molecular line emission of embedded, accreting 
massive protostars. We focus on molecular lines because 
they are the only available means of following gas kine- 
matics, and therefore of establishing whether observed 
structures truly are disks. In § [21 we describe our compu- 
tational method for making this prediction. This section 
is rather technical, and is likely to be of greatest interest 
to computational specialists. In § [Sj we present simu- 
lated emission data for a selection of molecular lines. We 
also process the data through simple models to mimic 
the performance of next-generation telescopes. In §2] we 
consider the implications of our results for future obser- 
vations, including how these observations may be used 
to distinguish between competing models of massive star 
formation. We also discuss the limitations of our model 
and how these are likely to impact our predictions. Fi- 
nally, in §[5] we summarize our conclusions. 

2. COMPUTATIONAL METHODOLOGY 
2.1. Calculation of the Emergent Intensity 
2.1.1. Ray Tracing Through an Adaptive Grid 

We wish to compute the radiative intensity emerging 
from a region within which information on the density, 
velocity, and temperature fields is recorded on an adap- 
tive grid consisting of a series of levels I = 0,1, ... L on 
which the cell size is Aa;;. The coarsest, lowest resolution 
level is Z = 0, and the highest resolution level isl = L. To 
compute the intensity from this data, we first construct a 
grid of rays, representing different lines of sight from an 
observer at infinity, passing through the computational 
domain. The rays are evenly spaced in a rectangular grid 
on the plane at infinity, with adjacent rays separated by 
a distance Axl, where Axl is the length of the finest cell 
in the adaptive grid. Along each ray, we construct a list 
of the cells through which it passes, ordered so that cell 
is closest to the observer. At every point in space, we 
record in the list only the finest resolution cell covering 
that point. Figure [T] illustrates this procedure. 

From the cell list for a given ray, we construct the list of 
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Fig. 1. — Illustration of rays (dashed lines) passing tlirougli an 
AMR grid. Cell edges are indicated by thin solid lines, and grid 
boundaries are indicated by thick solid lines. Cells are numbered 
by {column, row) from the lower left corner of a grid, starting with 
cell 0, so the point marked by the asterisk is in cell (1,3) of level 
1, grid 0, and also cell (1,2) of level 0, grid 0. The cell list for 
ray A is level 0, grid 0, cells (0, 0), ... , (6, 0). The list for ray B 
is level 0, grid 0, cell (0, 1); level 1, grid 0, cells (0, 1), . . . , (3, 1); 
level 0, grid 0, cell (3, 1); level 1, grid 2, cells (0, 1), . . . , (3, 1); level 
1, grid 0, cell (6, 1). The cell list for ray C is level 0, grid 0, cells 
(0,4), (1,4); level 1, grid 1, cells (0, 0), . . . , (3, 0); level 1, grid 2, 
cells (0, 6), . . . , (3, 6); and level 0, grid 0, cell (6, 4). 

densities po, . . . pn, line-of-sight velocities vg, . . . ,Vn, and 
temperatures Tq, . . . ,T„ through which the ray passes, 
with subscripts increasing away from the observer, and 
velocities defined so that positive corresponds to motion 
away from the observer. We also have the list of the 
lengths of the ray segments dsg , . ■ ■ dSn intersecting each 
cell. We treat the density, velocity, and temperature 
fields as step functions, so that at distance s along the 
ray such that Si < s < s^+i, the density, velocity, and 

temperature are pi, Vi, and Ti, where Si = J2]=o dsj, and 
s = corresponds to the point at which the ray enters 
the computational domain. 

2.1.2. Emissivities and Absorption Coefficients 

Now that we know the density p, line-of-sight veloc- 
ity V, and temperature T along each ray, we must use 
these to compute the emissivity and absorption coeffi- 
cient along the ray. We consider two emission and ab- 
sorption mechanisms. The first is dust thermal emis- 
sion a nd absorption. Following IChakrabarti fc McKeel 
(12005), we adopt a modified version of the 'Drain^ 
( 20()3a,.b>c,) opacity law (w-hich i s an updated version of 
the IWeingartner fc Drain"il2001l law). We set the dust 
specific absorption coefficient k'^^^^ at frequency v equal 
to that predicted by Draine's model for Ry = 5.5 in gas 
with a temperature above 110 K. In colder gas we set 
the absorption coefRcient to double that of the Draine 
model to account for the extra opacity provided by ice 
mantles on dust grains. Note that k'^™^ is the opacity per 
unit gas mass, not the opacity per unit dust mass. Thus, 
the total absorption coefficient is simply a^"*'* — n'^^^^p, 
and since the dust emits thermally, the emissivity is 
j^""* = a;^""'B^(T), where B^{T) is the Planck func- 
tion. At the radio and submillimeter wavelengths with 
which we are concerned, scattering is negligible and we 
need not distinguish between absorption and extinction. 

The second source of emission and absorption is molec- 
ular lines. Consider a line emitted by a species X transi- 
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tion from an upper state u to a lower state I. The line is 
at a frequency i^ui in the frame comoving with an emit- 
ting atom. In this paper we limit ourselves to species 
and levels that are in local thermodynamic equilibrium 
(LTE). Since the density in the massive cores we wish 
to examine is everywhere larger than ~ 10^ cm~^, and 
in the disk it is typically ~ 10^*^ cm~'^, this is not a sig- 
nificant restriction. We discuss the validity of the LTE 
approximation for individual species in § [31 For a species 
in LTE, the number density of molecules of species X in 
state I is 



ni = fxgi — 

fJ-H 



exp 



El 



Q{T) 



(1) 



where fx is the abundance of species X relative to hy- 
drogen, gi is the statistical weight of state is the 
mean gas mass per H nucleus (= 2.34 x lO"^"* g for a gas 
of hydrogen and helium in the standard cosmic abun- 
dance) , El is the energy of state I relative to the ground 
state, and Q{T) is the partition function for species X. 
The corresponding absorption coefficient at frequency v 
is 
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where i3/„ is the Einstein coefficient for absorption from 
state I to state m, E^i is the energy difference between 
states u and I, and (p{v) is a function describing the line 
shape. In LTE, the source function is simply the Planck 
function, so the emissivity 
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To determine the line shape function, we assume that 
the gas has a Maxwellian velocity distribution, and that 
both the bulk velocity and the velocity dispersion are 
piecewise-constant along the line of sight. We ignore 
the intrinsic width of the line relative to the thermal 
and bulk-motion broadening. If the bulk velocity of the 
gas is V, then a photon of frequency v in the observer's 
frame has frequency v' = (1 + (3)v in the frame comoving 
with the gas, where (3 = v/c. The velocity dispersion 
of species X along the line of sight is = ^JkBT/mx, 
where mx is the mass of particles of species X. The cor- 
responding dispersion in frequencies is = {(Jy/c)uui. 
Combining the bulk velocity and the velocity dispersion 
of the gas, the line shape function is 



: exp ■ 



2al 



(4) 



In practical calculations it is generally more convenient 
to work in velocity offsets relative to the mean velocity of 
the gas than in frequencies. If we define an observation 
at velocity Vohs to mean an observation at a frequency 
V = (\ — Vohs/cji'ui, then we may alternately write the 
line shape function as 



1 



exp 



{V - Wobs)^ 



2a?, 



(5) 



For all the calculations presented in this paper, we per- 
form our computations over velocities ranging from —20 
to 20 km s^^ relative the velocity of the primary star (see 



i) l2.3p . with a velocity resolution of 0.1 km s^^. Since 0.1 
km s~^ smaller than the thermal velocity dispersion of 
the molecules we discuss below even in the coolest parts 
of our computational domain, this resolution is sufficient 
to sample the emission well. 

Our assumption that the velocity is piecewise constant 
is potentially problematic, since in reality there are ve- 
locity gradients along the line of sight that correspond 
to velocity differences from one side of a cell to another 
that are not necessarily small compared to the thermal 
sound speed of the gas. Failure to account for these gra- 
dients in our calculations causes the line profile we com- 
pute emerging from a given pixel to be somewhat jagged. 
However, we are interested in aggregate quantities that 
average over many cells and not in detailed line shapes. 
Since the jaggedness in our spectra does not significantly 
change integrated quantities such as the total intensity 
or intensity-weighted mean velocity at a given pixel, we 
do not attempt to correct for it by including a treatment 
of line of sight gradients. 

2.1.3. Solution of the Transfer Equation 

Now that we have determined the emissivity and the 
absorption coefficient, the final step is to solve the trans- 
fer equation, 

-T- ^ -jiy + a^L, (6) 
as 

where is the radiative intensity at frequency i/ trav- 
eling in the —s direction, and ju = jjj"*'* -I- jj,™'' and 
ai, — aj^"*** + are the total emissivity and absorp- 
tion coefficient from dust and lines. Note that taking I^, 
to be positive in the — s direction inverts the usual sign 
convention, and this is the reason our transfer equation 
also has the reverse of the usual sign convention. 

Fortunately, since the emissivity and the absorption 
coefficient are functions simply of density, velocity, and 
temperature, and these are piecewise-constant on our 
grid, ji, and ai, are likewise piecewise-constant. This 
makes solution of the transfer equation trivial. For ray 
segment i, within which the emissivity is j^^i and the ab- 
sorption coefficient is a,y^i, integrating the transfer equa- 
tion ^ from Si (the boundary between segments i and 
i -I- 1) to Si_i (the boundary between segments i — 1 and 
i) gives 

Ii^isi-i) = exp{~T^^i)I^{si) + [1- exp{-T^^i)]B^^i, (7) 

where r^^i = au,idsi is the optical depth through segment 
i, B^^i is the Planck function evaluated at the tempera- 
ture in segment i, and Iv{si) is the intensity entering the 
segment. Combining a series of integrations of the form 
([7]), and taking the intensity entering the computational 
domain from the direction opposite the observer to be 
zero, the intensity emerging to reach the observer from 
s = is simply 

^.(0) = 51 '^^P ^ H V« [1 - exp(-v,)] B,,,. (8) 

2.2. Input Molecular Data 

We obtain the molecular data required for the compu- 
tation described in § 12.11 from the Jet Pr opulsion Lab- 
oratory Molecular Spectroscopy Catalog^ (jPickett et al.l 

http://spec.jpl.nasa.gov/ 
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I1998D . For each species, the catalog provides the value 
of the partition function at temperatures of 9.375, 18.75, 
37.5, 75, 150, 225, and 300 K. We compute our partition 
function at other temperatures by fitting a piecewise- 
powerlaw function to the catalog values. At tempera- 
tures above 300 K, we extrapolate using the same pow- 
erlaw index as between 225 and 300 K. We take the abun- 
dances of each species we use from values available in the 
literature, as described below in § [3) 

2.3. Input Simulation Data 

We apply our radiative transfer calculations to the den- 
sity, velocity, and gas temperature fields produced by the 
simulations of KKM07. Wc only briefly review the sim- 
ulation methodology and initial conditions here, and re- 
fer re aders to KKM07 and th e companion methodology 
paper iKrumholz et al.l (|2007a[ ) for a full description. Al- 
though the radiative transfer method we describe in this 
paper is applicable to any grid-based simulation, all the 
predicted emission line data we present here comes from 
applying this method to run lOOA of KKM07, so we focus 
on it in our description. This run begins from a spherical 
core with a mass of 100 Mq and an initial radius of 0.1 
pc. The core is initially centrally condensed, with den- 
sity profile p~^'^ at large radii, and a maximum initial 
density of lO"^'^ g cm"^ at radii < 38.4 AU. The core 
has an initial turbulent velocity field with a velocity dis- 
persion of 1.7 km s""'^ chosen to put the initial core into 
approximate hydrostatic balance. The turbulent power 
spectrum is P{k) oc k~^. The initial temperature in the 
core is 20 K, and the radiation boundary condition in the 
simulation is that there is a uniform radiative flux equal 
to that of a 20 K blackbody radiation field entering the 
computational domain. Radiation generated within the 
domain may escape freely at the domain edge. 

KKM07 evolve the core using the Orion code, which 
solves the equations of radiation-hydrodynamics on an 
adaptive mesh in the gray, flux-limited diffusion approx- 
imation. The simulation domain is a cube 0.6 pc on side, 
and the peak resolution is 7.5 AU, a dynamic range of 
16, 384 cells per linear dimension. KKM07 report the re- 
sults of the first 20 kyr of evolution, but we have contin- 
ued the simulations using the same code, and we present 
simulated maps based on the state of the core after 27 
kyr of evolution. At this time, the core has formed a 
primary protostar of 8.3 Mq, surrounded by a disk 3 — 5 
Mq in mass, and 500 — 1000 AU in radius. The ex- 
act mass and size of the disk are somewhat indetermi- 
nate, since in the simulation the disk has no sharply- 
defined edge. The central star has completely burned 
its accumulated deuterium reserve but has not yet con- 
tracted onto the main sequence. Through a combination 
of Kelvin-Helmholtz contraction and accretion power it is 
producing a time-averaged luminosity of approximately 
XO^-S _ lo^ Tjq. xhe luminosity and mass of the central 
object make the simul ated core an analog of systems such 
as IR AS 20126-f4104 (jZhang et al][l999l : ICesaroni et all 
|2005( ). In addition to the primary star, there are two 
other stars in the simulation domain that are 1.8 and 
0.05 Mq in mass, neither of which is producing much lu- 
minosity in comparison to the primary star. We show the 
column density in the simulation at this time in Figure 
[3 Note that, viewed in either orientation, the column 
density through the disk reaches almost 1000 g cm~^. 



Before using the simulation data for our radiative 
transfer calculations, we modify it in two ways. First, 
in the simulation KMM07a impose a constant pressure 
boundary condition by surrounding the core with an am- 
bient medium within which the temperature is 100 times 
higher and the density is 100 times lower than the gas 
at the edge of the core. To avoid contaminating the pre- 
dicted intensities with contributions from this medium, 
we set the density to zero in any cell within which the 
density is less than 5 times the initial ambient medium 
density. Since this is a factor of 20 lower than the initial 
core edge density, this threshold should not eliminate a 
significant amount of core material. 

The second modification is to the temperature field. In 
any given time step there are usually a handful of cells 
where the iterative radiation solver in Orion has not con- 
verged, and the temperature is unphysical. While these 
cells have no significant effect on the dynamics, since 
they disappear after a single time step and never con- 
tain a significant amount of mass, unphysically warm 
cells can affect the predicted intensities. A single cell 
that is twice as hot as its neighbors can outshine them 
even if it contains very little mass, particularly in veloc- 
ity channels close to that cell's velocity. We therefore 
remove these unphyiscally warm cells as follows. The 
signature of the non-convergence is a rapid oscillation in 
temperature over a few cells. We detect this signature by 
flagging as suspect any cell where, along the line-of-sight 
for which we are computing the transfer, there is both a 
local maximum and a local minimum in the temperature 
within 6 computational cells. We replace the tempera- 
ture in such flagged cells with a linear interpolation of 
the temperatures of the neighboring non- flagged cells. 
Visual inspection of the modified temperature fields that 
we produce using this procedure shows that they agree 
well with what one would identify and correct by eye. 

As a final point, although in this paper we only present 
radiative transfer calculations for run lOOA of KKM07, 
we have processed other simulations from KKM07 us- 
ing the same techniques. The details of the predicted 
line emission obviously differ from run to run, since the 
underlying density and temperature distributions differ 
in detail. However, we find that there is no significant 
qualitative difference between the line emission predicted 
for run lOOA and that predicted for other runs at times 
when the primary star's mass is comparable to 8.3 Mq, 
the mass of the primary in rmi lOOA at the time we con- 
sider. 

3. PREDICTED MOLECULAR LINE EMISSION 

3.1. Choice of Lines 

In principle we can use the method described in § [2] 
to compute the predicted emission from any molecule in 
LTE. At the densities of ~ 10^ cm^'^ or more typical of 
the structures we are interested in tracing, this includes 
most species. However, selection of appropriate species is 
not trivial. It is only possible to obtain kinematic infor- 
mation if the emission at a given frequency is dominated 
by molecular emission rather than dust emission, or the 
two are at least comparable. Since the intensity emitted 
by the disk along a given line of sight cannot exceed the 
value of a blackbody, and if the dust is optically thick it 
will radiate as a blackbody, it is only possible for molec- 
ular emission to dominate if the dust is optically thin. 
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Fig. 2. — Column density in regions 10000 AU and 1000 AU in size (upper and lower rows) in two projections [left and right 
columns). Stars are shown in the 1000 AU images by plus signs. The 8.3 Mq star is at the origin of both plots. In the image 
on the left, the 1.8 Mq star is to the one to the left of the origin; in the right image, the 0.05 Mq star is the one on the far left. 
The overdense blob attached to the disk spiral arm, shown near the bottom of the bottom left image and the left of the bottom 
right image, is in the process of collapsing to form an additional fragment, but has not collapsed yet. Note that the drop in column 
density in the immediate vicinity of the 8.3 and 1.8 Mq stars is an artifact of the finite resolution of the simulation, not a true physical effect. 



This is a significant challenge: at 100 GHz, the opacity 
predicted by our dust model is k^J"^* = 5.6 x lO"** cm^ 
g~^, and the opacity in cold regions where grains have 
ice mantles will be 1.1 x 10"'^ cm^ g~^. (Again, note that 
this is the opacity per unit gas mass, not the opacity per 
unit dust mass.) Figure [2] indicates that the column den- 
sity in the inner disk, within ^ 100 AU of the primary 
star, falls in the range 100 — 1000 g cm~^. This means 
that the inner disk has optical depth of at least a few 
tenths, usually more, in dust, and thus is it is not pos- 
sible to obtain velocity information for it at 100 GHz or 
higher frequencies. 

This calculation suggests that the best tool for tracing 
very small scales is radio observations at frequencies well 
below 100 GHz, where the dust contribution is smallest. 
Unfortunately, the brightness temperature sensitivity to 
line emission for current and future radio telescopes such 
as the VLA and EVLA is much poorer than the sensitiv- 
ity of current and future submillimeter telescopes such 
as PdBI, the SMA, and ALMA. Thus, the hot inner disk 
is likely to be the only part of the structure around our 



protostar detectable at radio frequencies using the high 
resolutions required to probe the structure of the disk, 
and this will likely require very long integration times. 
Observations at frequencies of 100 GHz or more are far 
more efficient for tracing somewhat lower column density 
structures on larger scales. 

Based on these considerations, we present simulated 
observations for the lines of two molecules. At radio 
frequencies we use the inversion transitions of ammo- 
nia (NH3), and in the submillimeter regime we use some 
of the strongest rotational transitions of methyl cyanide 
(CH3CN). Both molecules have previously been used to 
trace structures around massive stars on larger scales and 
at lower resolutions with the curre nt generation of tele- 
scopes (see the rece nt r eviews bv ICesaroni et al.l l2006l 
120071 iBeuthol 120061 and iBeuther et al.ll2007D r and thus 
are known to be present and excited in massive protostel- 
lar envelopes. We simulate observations for distances of 
0.5 kpc and 2 kpc. The former is roughly the distance to 
the Orion molecular cloud, the nearest region of massive 
star formation, while the latter is a more typical for ob- 
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Fig. 3. — For the disk seen face-on, velocity-integrated brightness temperature emitted in the indicated NH3 inversion transitions {upper 
row), and integrated brightness temperature detected in a simulated observation of these lines using the EVLA (lower row). In the 
simulated observation, we show as black any pixel for which there is no detection at the 3cr(= 144 K) level in any velocity channel. The 
beam size for the simulated observation is shown by the hatched region in the lower right-hand corner. Contours start at 500 K km s~^, 
and increase by factors of 2 per contour. 



servations of massive star-forming regions. We assume 
that all telescope beams are Gaussian, and do not in- 
clude the removal of large-scale power by interferometers 
in our models. 

3.2. Ammonia 

We first compute the radiation intensity emerging from 
a massive protostellar core in the {J,K) = (2,2), (4,4), 
(6, 6), and (8, 8) inversion transitions of NH3, at frequen- 
cies of 23.7226, 24.1394, 25.0560, and 26.5190 GHz. The 
upper levels for these transitions have temperatures of 
65.0, 201.1, 408.6, and 687.4 K above the ground state, 
so they probe a wide range of temperatures, including 
those that are expected to occur only in the inner disk 
very close to the star. At temperatures > 20 K, the crit- 
ical densities for the upper states of these transitions are 
all ^10^ cm~^. Since the density everywhere in our core 
is si 10^ cm^"^, our assumption of LTE for these levels is 
well-justified. (To compute the critical densities, we use 
the molecular data a vailable in the Leid en Atomic and 
Molecular Database^. ISchoier et al.ll2005l since the JPL 
database does not include the collision strengths required 
to compute critical densities.) 

The abundance of NH3 in dense molecular cores is 
fairly uncertain. In low mass, cold cores observa- 
tions show__thBt ammonia abundance rises with den- 
sity (|Tafalla et al.ll2004allbl ). reaching values of 10"^ at 
the high densities found in our core. Observations 
of massive star-forming regions generally find signifi- 
cantly lower abundances o f 10~^ — 10"'' (jHariu et al.l 
1X9931: Irieftrunk et al.lll998l ). However, these abundances 

^ http:/ /www. strw.lcidcnuniv.nl/~moldata/ 



are based on single-dish observations of distant regions 
that average over size scales of ^0.5 pc, and thus 
presumably include a great deal of low density mate- 
rial w ith correspondingl y low ammonia abundance. In- 
deed, iHariu et al.l |1993) find that their derived ammo- 
nia abundances decrease systematically with distance, a 
result they attribute to a decrease in the NH3 filling fac- 
tor as their fixed-angular size beam averages over larger 
and larger spatial scales. It therefore seems likely that 
the NH3 abundance in our disk is significantly higher 
than the values derived from large-scale averaging, and 
we adopt the larger value of 10~^ determined for low 
mass cores as our fiducial abundance. Fortunately, our 
results are fairly insensitive to this choice, because the 
NH3 line is optically thick through the bulk of the disk 
unless the abundance is more than an order of magnitude 
lower than our fiducial value. 

Finally, note that our calculations are only for the main 
part of the NH3 line, not the hyperfine satellites. These 
satellites are intrinsically weaker by a factor of ~ 14. 
While this would decrease their optical depths and thus 
make them more faithful tracers of kinematics, as we 
show below even the optically thick NH3 main line emis- 
sion is difficult to detect even in long integrations. These 
weaker lines would be harder still to observe. 

Figures [3] and [5] show the velocity- integrated intensity 
in the four NH3 lines for our core at a distance of 0.5 
kpc viewed nearly face-on and nearly edge-on, in the ori- 
entations shown in the left and right columns of Figure 
[21 In both Figures, the upper row shows the velocity- 
integrated brightness temperature as a function of posi- 
tion relative to the primary star. As one might expect, 
the low-temperature (2, 2) line comes from a large re- 
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Fig. 4. — For the disk seen edge-on, velocity-integrated brightness temperature. See the caption of Figure|3]for details of what is shown 
in each panel. The dashed lines in the simulated observations indicate the orientation of the disk, inclined 15° with respect to the x axis 
in the observation. 
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Fig. 5. — Brightness temperature in the indicated NH3 inversion transitions as a function of position and velocity, with position measured 
relative to the star along the dashed line indicated in Figure |4] We show both the true brightness temperature {upper row) and a simulated 
EVLA observation {lower row). The contours start at 50 K and increase by factors of 2. In the simulated observation, we report for any 
pixel in which the signal is smaller than the 3a{= 144 K) noise level. The angular and velocity resolution of the observations are indicated 
by the error bars in the lower right-hand image. The thick dashed lines indicate the velocity profile of a Keplerian orbit around a star of 
8.3 Mq, offset by a constant velocity of —3 km s~^ relative to the zero velocity in the image. 



8 



Krumholz, Klein, & McKee 



gion extending hundreds of AU away from the central 
star, while the (8, 8) line, which requires temperatures of 
nearly 700 K to be excited, comes almost entirely from 
a region within 100 AU of the primary and secondary 
stars. Also note that the disk around the 1.8 Mq sec- 
ondary appears quite distinctly, and that there is also a 
hot spot visible around the 0.05 Mq star. 

The lower rows of Figures [3] and |4] show simulated ob- 
servations of the disks using the EVLA phase I, which 
at the frequencies of the NH3 inversion lines is projected 
to have a FWHM resolution of 0."12 in the A configura- 
tion. We adopt a velocity resolution of 4 km s~^ for the 
simulated observation. The estimated 1 a thermal noise 
of the EVLA in A configuration for a 12 hour integration 
in K- band with this bandwidth is 48 K (jPerlev fc RupenI 
1200 It ). We take the la noise in the integrated intensity to 
be the noise in a channel times the channel width, 192 K 
km s~^. Note that, although the EVLA is capable of far 
finer velocity resolution, simulated observations with sig- 
nificantly smaller velocity channels produce images with 
very low signal to noise ratios; since 4 km s~^ resolu- 
tion is sufficient to map the rotation of the inner disk, 
we adopt it. With these observational parameters, the 
system is reasonably well resolved, and the presence of a 
flattened disk is obvious in the edge-on observation. The 
secondary is also clearly resolved as a separate peak in 
the in the edge-on images, although the material around 
it is not clearly shown as a disk due to beam smearing. 
The tertiary star is outside the box shown, and the gas 
around it is below the sensitivity limit of the observation 
in both orientations. 

Figure [5] shows the true brightness temperature and a 
simulated EVLA observation as a function of velocity and 
position along the disk viewed edge-on, indicated by the 
dashed line in FigurelH As in the total intensity map, the 
disk is reasonably well-detected and marginally resolved 
in all four transitions. The velocity profile of the disk is 
well-fit by a Keplerian orbit around an 8.3 Mq star, the 
true mass of the central object. Interestingly, however, 
in order to obtain a good fit to the true emission, the 
mean velocity of the orbit must be shifted by —3 km 
s""'^ relative to the velocity of the star, which is zero on 
the axes shown. However, the offset persists even in the 
true brightness temperature map. Its physical origin is in 
the fact that the star is really moving with respect to the 
disk, due to the large mass of the disk relative to the star. 
This is part of the gravitational instability mechanism 
for angular momentu m transport which is opera t ing in 
the disk, as discussed iLaughlin fc Bodenheimeil ()1994D 
and KKM07. In the simulated observations a somewhat 
larger offset would be required to obtain a good fit, but 
the difference is below the velocity resolution limit of the 
simulated observation. 

Whether this offset between the central velocity of the 
disk and the stellar velocity is detectable observationally 
depends on whether it is possible measure the velocity of 
the star. Since the star is deeply embedded, one cannot 
measure its velocity directly. However, if such an em- 
bedded star produces a relatively small ultraviolet fiux, 
it may generate a hypercompact HII region which will 
be detectable in radio emission through the intervening 
gas. During the hypercompact phase when ionized gas 
is confined very close to the stellar surface, the ionized 
gas velocity should presumably be close to the stellar 



velocity. The 8.3 Mq star we model here probably pro- 
duces too little ionizing fiux to generate a detectable HII 
region, however. 

We also compute simulated observations for a distance 
of 2 kpc. However, at this distance the entire region of 
NH3 emission detectable in a 12 hour integration falls 
within a single beam pointing at the EVLA's highest 
resolution. For this reason, we do not discuss these sim- 
ulated observations further. 

3.3. Methyl Cyanide 

We next compute the emission from our disk in 
four rotational transitions of CH3CN: (12,0,0,13) 
(11,0,0,12), (12,3,0,13) ^ (11,3,0,12), (12,5,0,13) ^ 
(11,5,0,12), and (12,0,2,13) ^ (11,0,2,12). Here the 
quantum numbers are given in the JPL catalog format 
(N^ K,v, J) for symmetric rotors with vibration, where 
N is the total angular momentum of the molecule ex- 
cluding electron spin, K is the projection of this onto the 
molecule's symmetry axis, v is the vibrational quantum 
number, and J is the total angular momentum includ- 
ing electron spin. These transitions occur at frequen- 
cies of 220.7472, 220.7091, 220.6411, and 221.3942 GHz, 
and their upper states lie 69.0, 133.3, 247.6, and 594.6 
K above the ground state. We choose these transitions 
because they are among the strongest for the molecule, 
with Einstein A's of almost 10^'^ s^^. 

Unfortunately collision strengths for the CH3CN levels 
we examine are not available in the astrophysical litera- 
ture, so we cannot compute critical densitie s. However, 
collision strengths for lower CH3CN levels ()Pei fc Zengj 

IT995I) give critical densities ^ 10^ cm ^ even for the 
strongest radiative transitions, as do collision strengths 
for molecules of similar complexity such as methanol 
(CH3OH). (Data for methanol a re from the Leide n 
Atomic and Molecular Database, ISchoier et al.l l2005l ) 
Since we are concerned with emission from the disk and 
structures around it, where typical densities are gener- 
ally 10^ cm~^ or more, it seems reasonable to assume 
that methyl cyanide molecules are in LTE. 

The abundance of CH3CN is, unfortunately, even less 
well understood than that of NH3. The formation of 
methyl cyanide is driven by complex organic chemistry 
that occurs at the comparatively high temperatures of 
hot cores like the one in our simulation. There is strong 
evidence that the abundance of nitrogeno us organic 
molecules is time-dependent (jBeuthed l2006l ). Rather 
than attempt to add a chemistry model to our radiative 
transfer calculation, we adopt a fiducial CH3CN abun- 
dance of 10^* independent of density. This is roughly 
cons istent with the values measured by iHatchell et al.l 
(|199^ along lines of sight toward 14 ultracompact HII 
regions. We discuss how variations in the abundance will 
affect our results below. 

Figure [5] shows the velocity- integrated brightness tem- 
perature in the four CH3CN transitions seen edge-on. 
As with Figures [3] and IH the upper row shows the true 
brightness temperature. Since dust thermal emission 
dominates essentially everywhere when integrated over 
the 40 km s^^ interval for which we compute the output, 
this is effectively a dust continuum image integrated over 
a 40 km s^^ channel. To isolate the molecular line emis- 
sion (and in some cases absorption) we subtract from 
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Fig. 6. — Velocity-integrated brightness temperature emitted in the indicated CH3CN rotational transitions over a 40 km s~^ velocity 
range centered about the center of each line. We show the true total emission {upper row), the true continuum-subtracted emission {second 
row from the top), the total emission as observed with ALMA {third row from the top) at 0.5 kpc, the continuum-subtracted emission as 
observed with ALMA {fourth row from the top), and simulated ALMA observations of the total and continuum-subtracted emission at 2 
kpc {fifth and sixth rows, respectively). The upper level of the transition shown in each image is indicated in each column. The size of the 
ALMA beam at 0.5 kpc and 2 kpc is indicated in the lower right corners of the fourth and sixth rows. In the true continuum-subtracted 
emission, a black pixel indicates a position where the line is seen in absorption rather than emission. In the simulated observations, black 
pixels correspond to positions where there is no detection at the 3cr(= 33 K) level in any velocity channel. For the total emission images 
{odd rows), the contours start at 50 K km s~^; in the continuum-subtracted images {even rows), they start at 25 K km s"-"^. In all cases 
they increase by factors of 4 per contour. Note the difference in scales from Figures [3] and |4l and also note that the color scales are different 
for the total and continuum-subtracted images. 
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each bin the emitted intensity in the —20 km s^^ bin, 
which contains a neghgible contribution from Une emis- 
sion, in order to produce a continuum-subtracted image. 
We show this in the second row of Figure [6l The inner 
disks around each star are mostly seen in absorption or 
weak emission, while the larger-scale extended structure 
shows strong line emission. 

Note that our continuum-subtraction procedure is only 
a rough approximation to the continuum subtraction pro- 
cedure used in real interferometric observations, which 
involves fitting the observed visibilities in the uv-domain 
with polynomials. However, the qualitative conclusions 
we draw about the effects of continuum-subtraction do 
not depend in detail on this point. 

The bottom four rows show simulated observations of 
the total and continuum-subtracted images made with 
ALMA, at distances both 0.5 kpc and 2 kpc. For the 
simulated ALMA observations, we adopt a velocity res- 
olution of 0.5 km s^^, an angular resolution of 0."1, and 
assume 1000 seconds per pointing. With these values, 
the la noise level is roughly 11 K per channel for each 
line. (We compute the noise level using the ALMA sen- 
sitivity calculator'* for 50 antennas.) 

As expected, in all four transitions the inner disk, 
within ^ 100 AU of the primary star, is detected in ab- 
sorption of the dust continuum by molecules in cooler 
foreground gas. Absorption information is difhcult to 
use to determine disk kinematics, because there is no 
easy way to separate absorption arising from a cool sur- 
face layer of the disk, which would trace its motion, from 
absorption arising from foreground gas elsewhere in the 
core that is not part of the disk. As a result, we will not 
attempt to use the line detected in absorption to trace 
kinematics of the inner disk. On the other hand, the 
outer disk at larger radii is well-detected and resolved 
in emission in two of the four lines, and marginally de- 
tected in the third. The highest temperature line pro- 
vides no significant signal because the region of the core 
hot enough to emit it is largely confined to the region 
where the dust is optically thick and the line is detected 
primarily in absorption. 

This calculation also illustrates the importance of res- 
olution and distance. With infinite resolution, as indi- 
cated in the second row of Figure [6l the entire disk is 
visible in line emission. With lower resolution, however, 
it is not possible to subtract off the dust emission faith- 
fully and recover the line emission on size scales compa- 
rable to the beam size or smaller. The optically thick 
continuum emission of the inner disk and the gas imme- 
diately around the star raises the emitted intensity to the 
blackbody intensity, and since this is the maximum pos- 
sible emitted intensity in LTE, the line cannot be seen 
above it. Consequently, less and less of the disk becomes 
observable as the distance, and thus the physical scale 
probed by a beam of a given angular size, increases. In 
order to probe the inner disk as close to the star as pos- 
sible, one therefore wants the highest possible angular 
resolution. 

In the transitions where the disk is detected, the strong 
spiral structure of the disk is readily apparent at either 
distance. One can confirm that the structure is a spiral 

http:/ /www. eso.org/projects/alma/science/bin/ 
sensitivity, html 



arm in a disk rotating about the central object by ex- 
amining its velocity structure. Figure [7] shows the true 
value and a simulated ALMA observation of the centroid 
velocity of the continuum-subtracted line emission. We 
define the centroid velocity as the continuum subtracted- 
brightness temperature- weighted velocity average at each 
point. In computing the average for the true brightness 
temperature, we only include pixels for which the contin- 
uum subtracted-brightness temperature is positive, for 
the true value map, or is greater than the 3cr noise level 
of the beam, for the simulated observation. 

Although it is not directly possible to compare the 
rotation profile shown by Figure [7] with the Keplerian 
speed, since the disk is clearly non-circular and the dis- 
tance from the central object to the point in the arm 
from which emission arises is not clearly determined, the 
simulated observation clearly does reveal rotation con- 
sistent with what one would expect for a spiral arm in a 
disk. At 500 - 1500 AU from the primary 8.3 Mg star, 
the Keplerian speed is 2.2 — 3.8 km s~*, consistent with 
the observed velocities for a distance of 0.5 kpc. At an 
assumed distance of 2 kpc, the velocities are considerably 
lower, due to blending of material at different velocities. 
Also note that both the true brightness temperature map 
and the simulated observation reveal significantly more 
negative velocity than positive velocity gas; the entire 
structure is shifted away from zero velocity. This is the 
same offset seen in the NH3 observation on smaller scales. 
However, it is important to remember that the centroid 
velocities are not necessarily perfectly reliable tracers of 
the true velocity, particularly toward the denser parts of 
the spiral arm, because the lines becomes optically thick 
in places. 

Finally, we return to the question of how much our 
results depend on the assumed abundance of methyl 
cyanide. To help explore this question, we calculate emis- 
sion from the same four lines, adopting abundances a fac- 
tor of 10 higher and lower than our "best guess" value of 
10~*. For simplicity, we only compare calculations at 0.5 
kpc, not 2 kpc. Not surprisingly, we find that our results 
depend fairly strongly on the abundance. With an abun- 
dance of 10^^, we find almost no detectable continuum- 
subtracted line emission using the same sensitivity and 
resolution as for our calculations with an abundance of 
10~*. In contrast, figures [5] and [5] show the velocity- 
integrated intensity and centroid velocity as a function of 
position for an abundance of 10~^. In this case, we find 
line emission within ALMA's sensitivity arising from a 
significantly larger region, and find detectable emission 
after subtracting the continuum from all four lines, al- 
though the detection for the highest temperature line 
is over a fairly small area. However, the emission in 
this case traces the kinematics less faithfully than for an 
abundance of 10~^, particularly for the lower levels, be- 
cause the line is quite optically thick in many regions. For 
example, the peak positive velocity in the (12,0,0,13) 
and (12,3,0,13) lines is noticeably lower than the peak 
positive velocity in those lines for an abundance of 10^^, 
or than the peak positive velocity in the (12, 5, 0, 13) line 
for an abundance of 10~^. 

4. DISCUSSION 
4.1. Implications for Future Observations 
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Fig. 7. — Centroid velocity as a function of position in the indicated CH3CN rotational transitions, after subtracting off dust continuum 
emission. We show the true centroid velocity {top row) and the result of a simulated observation with ALMA assuming a distance of 
0.5 kpc (middle row) and 2.0 kpc (bottom row). The contours are separated by intervals of 2 km s~^. The thick contours correspond to 
negative velocities, and the thin contours correspond to positive velocities. 



Our analysis shows that the ~ 1000 AU-sized disks 
around young, accreting massive protostars should be 
detectable and at least marginally resolvable with the 
EVLA and ALMA at distances of 0.5 kpc, the distance 
to the nearest region of massive star formation. Perhaps 
most excitingly, these detections should unambiguously 
reveal rotation at velocities comparable to the Keplerian 
speed, and should also reveal the presence of large spiral 
arms, or any other large-scale non-axisymmetry. They 
will therefore enable us for the first time to probe the 
structure of the accretion disks that feed gas onto massive 
protostars. However, these observations will be costly. 
On the ~ 0." 1 size scales resolvable by these telescopes, 
the brightness temperature is hundreds of Kelvin, and 
obtaining strong detections of these signals will require 
tens of hours of integration. 

Our calculations also reveal some of the strengths and 
weaknesses of radio versus submillimeter observations. 
The column density in the inner parts of disks around 
massive protostars can reach of order 1000 g cm^^, 
viewed either edge-on or face-on. This column density 
is high enough to make dust optically thick at submil- 
limeter wavelengths, making it impossible to detect line 
emission and trace disk kinematics for the inner parts of 
disks using submillimeter observations. Only radio ob- 



servations are able to trace the kinematics in the inner 
parts of disk. However, because existing and planned 
radio telescopes have relatively little sensitivity to lines, 
they are only likely to be able to detect the very inner 
parts of disks with reasonable integration times. In con- 
trast, the greater sensitivity of submillimeter telescopes 
makes them excellent tools for tracing the outer, lower 
column density parts of disks. However, one does require 
high angular resolution in order to minimize the amount 
of the outer disk concealed by the bright, optically thick 
dust emission from the inner disk. 

Our results on the optical thickness of massive disks 
are in interesting contrast with observations of disks 
around low mass protostars. Models of these disks based 
on multi-wavelength observations generally indicate that 
dust reaches optical depths of unity only at frequencies 
> 350 GHz and radii of < 50 AU fe.g. lLav et allflOfll . 
The difference in our results arises simply from the much 
greater column densities that massive disks can reach. 
Even class low mass protostellar disks typically have 
S ^ 100 g cm~^ at radii larger than 100 AU, an order 
of magnitude lower than the column densities reached in 
our disks at similar radii. 

4.2. Distinguishing Models of Massive Star Formation 
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Fig. 8.— 

One compelling reason to search for disks around mas- 
sive protostars, and to explore their properties, is to 
distinguish between competing models of massive star 
formation. The tw o dominant models are core accre- 
tion (Mc Kee fc Tan 2002, 2003, KKM07), in which mas- 
sive stars form by the collapse of massive gas cores di- 
rectl y to massive stars, and collision / competitive ac cre- 
tion (jBonnell et al.lll998L[200llBonnell fc Batell2005( l. in 
which all stars are born with masses of ~ 0.5 Mq, and 
massive stars reach their final mass by a combination of 
runaway gas accretion and star-star collisions in the cen- 
ters of stellar clusters. The evidence for each of these 
models have re cently been review e d by a number of au - 
thorsj including'KrumhqU ()2006a[ ) , iBonnell et all (|2007[ ) , 
andlBcuther ct al. (200^^ 

High resolution observations of the environments of 
massive protostars offer a strong method of distinguish- 



ing between the models, and more generally a method of 
determining the nature of the massive star formation pro- 
cess. In particular, the coalescence/competitive accre- 
tion model predicts that essentially all stars larger than 
3 Mq should suffer close encounters with other stars, 
truncating their dis ks to sizes smaller than ^ 30 AU 
(jBonnell et al.l [20031 ). The core accretion model, as ex- 
emplified by the KKM07 simulation, predicts that mas- 
sive stars should have disks a few hundred to a few thou- 
sand AU in size, with masses that are typically tens of 
percent of the mass of the primary star. These disks 
may also have other pr otostars embedded within them or 
around them (see also lKratter fc Matzneil [20061 ). which 
may or may not be detectable as "hot spots" in the disks, 
depending on the separation and orientation. Thus, the 
detection of structures such as those presented in this 
paper would be very strong evidence for core accretion. 



Line Emission From Massive Disks 



13 



-7.0 



-4,0 



-i,a 



2.D 



=1 



, , {^"^ , , 

-1.S0-D.7 D.75 -1. 50-0.75^ 8.0^^ 0.75 -1 .50-D.7 D.75 - 1 .50-0.75^ 8.0&^ 0.75 1 .51] 




Fig. 9. — Same as the top four rows of Figure[7l but with a CH3CN abundance of 10 ^. 



while non-detections with ALMA and the EVLA would 
be strong evidence for competitive accretion. 

4.3. Limitations of Our Calculations 

Our predictions are subject to several limitations. 
Most significantly is our poor knowledge of the abun- 
dances of tracer molecules, and how these change in time 
and space. As discussed above, there is considerable ev- 
idence that the abundances of both NH3 and CH3CN 
vary depending on the ambient density and tempera- 
ture, and on the time history of density and tempera- 
ture e xperienced by a particular fluid element. iBeutheil 
()2006l ) describes recent observational work toward devel- 
oping a chemical evolutionary sequence for hot molec- 
ular cores, but these efforts are still preliminary. Con- 
sequently, our abundance estimates could easily be off 
significantly. This would have strong effects on observa- 
tions, as is clear from comparison of Figures [6] with [3 
and [7] with [5] Due to this limitation it is probably best 
to regard our predictions as examples of what lines in 
the radio and submillimeter regimes may look like. For 
a given source, the particular tracer molecules we have 
used may be much more or less abundant, in which case 
one should be able to obtain data comparable to what 
we present using alternative molecules with transitions at 
similar frequencies and upper level temperatures. In the 
absence of a reliable, predictive chemodynamical model, 
the only option may be to try many molecules until a 
useful one is detected. 

Another limitation is in the KKM07 simulation that 
we use to produce our predicted line emission. The de- 
tails of the simulation and its approximations and lim- 
itations are discussed in detail in KKM07, but two in 
particular are worth mentioning here because of their 
potentially large impact on our predicted line emission. 
First, the KKM07 simulations do not contain protostel- 
lar outflows, which are often the most prominent features 



in some tracers. These are obviously missing from our 
predictions, and thus we cannot comment on potential 
confusions in determining whether a given observational 
feature is tracing disks or outflows. 

Second, the KKM07 simulations begin with isolated, 
spherical cores, and thus are likely to differ from real 
cores on scales comparable to the initial core radius, 0.1 
pc. On these scales, the geometry of the core at for- 
mation and interactions between the core and its en- 
vironment are likely to have significant effects on both 
morphology and kinematics. For this reason, we do not 
expect our cores to reproduce observations on scales of 
that are significant fractions of 0.1 pc. However, while 
gas outside the core may affect the morphology and kine- 
matics of the core gas, it is important to point out that 
this gas will not in itself present a signiflcant source of 
emission or opacity. Massive cores are generally embed- 
ded in mol ecular clumps wit h column densities ~ 1 g 
cm^^ (e.g. iPlume et al.lll997b . so we expect this to be 
the typical column density of any intervening material 
between the core and the observer. At 100 GHz, this 
corresponds to a dust optical depth ~ 10"'^, completely 
negligible in either absorption or emission. The line opac- 
ity will also be negligible in the lines we examine in this 
paper, because the lower states for these lines have exci- 
tation temperatures ^ 60 K. Since the gas in molecular 
clum ps is typically at much lower temperatures of 20 
K (Pl ume et anil997t ). there wiU be a negligible number 
of molecules excited into the lower states of our chosen 
transitions, and thus almost no molecules capable of ab- 
sorbing photons coming from the core. 

5. CONCLUSIONS 

Using detailed radiation-hydrodynamic simulations of 
massive star formation as a starting point, we predict 
the molecular line emission of the disks around massive, 
embedded protostars. These disks are distinct from the 
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mass ive, rotating toroids discussed by (jCesaroni et al.1 
I2OO6I ) and others, in that they are Keplerian rather than 
sub-Keplerian, they are geometrically thin, and they are 
much smaller in both spatial extent and mass. We find 
that the disks emit line radiation at brightness temper- 
atures ranging from hundreds of kelvins within a few 
hundred AU of the central star, down to many tens of 
kelvins at distances of ~ 1000 AU. This emission will oc- 
cur over velocity channels 10 — 20 km s~^ wide. These 
brightness temperatures and bandwidths should make 
line emission from massive protostellar disks detectable 
and at least marginally resolvable with next generation 
telescopes such as ALMA and the EVLA. 

Disks around young, massive protostars have ex- 
tremely high column densities in their centers, and the 
resulting high dust columns require radio observations 
to probe structures and kinematics within 100 — 200 AU 
of the central star. At submillimeter wavelengths it is 
possible to image disks in dust thermal emission, but 
not to obtain kinematic information from lines. Line 
observations in the radio over distances of a few hun- 
dred AU should reveal clear evidence for disks rotating 
at speeds comparable to the Keplerian speed around the 
central object. The zero velocity for the Keplerian ro- 
tation curve may be offset relative to the velocity of the 
central star by a few km s~^ due to relative motion of the 
disk and star induced by strong gravitational instability 
in the star-disk system, although it will only be possible 
to detect this offset if one can observationally determine 
the stellar velocity. 

Submillimeter telescopes, with their greater sensitiv- 
ities, are more suited to mapping the outer, extended 
parts of disks between a few hundred and 1000 — 2000 
AU from the primary star. Observations on these scales 
should reveal if massive protostellar disks have signifi- 
cant spiral structure, as the simulations of KKM07 pre- 
dict they should. They should also show clear evidence 
for rotation. However, even in the outer disk, obtain- 
ing kinematic information requires observation of partic- 
ularly strong lines and/or abundant molecules in order 
for the emission to be detectable in reasonable integra- 
tion times. It will also require high spatial resolution so 
that as little line emission from the inner disk as possible 
is masked by bright continuum emission from close to the 
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protostar. 

Observations of disks around massive protostars can 
play a decisive role in deciding between competing mod- 
els for the formation of massive stars. The core accretion 
model predicts that massive protostars should be sur- 
rounded by disks ~ 500 — 1500 AU in size, with masses 
that are appreciable fractions of the mass of the cen- 
tral protostar. In contrast, the competitive competitive 
accretion / collision model predicts an absence of such 
massive, extended disks. Given the potential power of 
high-resolution radio and submillimeter observations for 
probing the nature of the star formation process, we en- 
courage all simulators and modelers of star formation to 
make predictions that can be directly compared to ob- 
servations using the next generation of telescopes. 
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